% dgdx_red = zeros(4,4);
% dgdx_red(:,1:2)=J(:,1:2);
% dgdx_red(:,3:4)=J(:,19:20);
% dfdx_red = [-1;-1;1;1];
% 
% lambda = dgdx_red'\(-dfdx_red);
% 
% dgdx_rem = J(:,3:18);
% 
% dfdx_rem = -dgdx_rem'*lambda;
% 
% dfdx = [dfdx_red(1:2);dfdx_rem;dfdx_red(3:4)];

u = bfu(size(bfu)*[0;1]).data;

u(3) = 0.002;

J = f_dgdu(u);

epsilon1 = 10^-3;

token = zeros(m*N,1);

for i=1:m*N
    if abs(u(i)) >= epsilon1
        token(i) = 1;
    end
end

dfdx_red = zeros(sum(token),1);
dgdx_red = zeros(n,sum(token));

count = 0;
for i=1:m*N
    if token(i) == 1
        count = count + 1;
        dfdx_red(count) = sign(u(i));
        dgdx_red(:,count) = J(:,i);
    end
end

lambda = dgdx_red' \ -dfdx_red; % uses pseudoinverse if not square
% attempts to minimize the reduced residual

% calculates full dfdx
if sum(token) > n % overdetermined case
    dfdx = zeros(m*N,1);
    count = 0;
    for i=1:m*N % -dgdx'*lambda does not recover dfdx_red if lambda is overdetermined
        if token(i) == 1
           count = count + 1;
           dfdx(i) = dfdx_red(count); % retrieves values from the reduced dfdx
       else
            dfdx(i) = -J(:,i)'*lambda; % creates new dfdx values from lambda and J
       end
    end
else % underdetermined case
    dfdx = -J'*lambda; 
end



residual = dfdx + J'*lambda; % calculates the full residual

residual_l2norm = norm(residual,2);
residual_l1norm = norm(residual,1);